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Abstract 

We study the synchronization of totalistic one dimensional cellular automata 
(CA). The CA with a non zero synchronization threshold exhibit complex non 
periodic space time patterns and conversely. This synchronization transition 
is related to directed percolation. We study also the maximum Lyapunov 
exponent for CA, defined in analogy with continuous dynamical systems as 
the exponential rate of expansion of the linear map induced by the evolution 
rule of CA, constructed with the aid of the Boolean derivatives. The synchro- 
nization threshold is strongly correlated to the maximum Lyapunov exponent 
and we propose approximate relations between these quantities. The value of 
this threshold can be used to parametrize the space time complexity of CA. 
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Cellular automata (CA) are discrete dynamical systems that may exhibit complex space 
time patterns It has been observed that CA may be synchronized by a stochastic 

coupling We study all totalistic CA with four, five and six neighbors. For the CA 

considered, we find that a synchronization threshold is reached critically and that all CA 
with complex non periodic space time patterns have a positive threshold and conversely. 
We also find a strong relationship between the synchronization threshold and the maximum 
Lyapunov exponent (MLE) of CA g|7[|. 

Let us start considering the following asymmetric coupling for a continuous one- 
dimensional map f(x) H 



with x = x(t), x' = x{t + 1) (idem for y) and < p < 1. The function f(x) depends in 
general by one (or more) parameter a; let assume that, for the chosen value of a, f(x) is 
chaotic with Lyapunov exponent A, and that x(0) ^ 2/(0). Then, x(t) is always different 
from y(t) for p = 0, while for p = 1 x and y synchronize in one time step. There exists 
a critical synchronization threshold p c for which both trajectories x(t) and y(t) become 
indistinguishable in the long time limit and 



In what follow we shall try to develop similar relations for CA. We begin with a brief 
review of the definition of maximum Lyapunov exponent for CA based on a linear expansion 
of the evolution rule. We then present a synchronization mechanism and show that the 
distance between two realizations goes to zero in a critical manner at p c . The numerical 
experiments show a relation between p c and the maximum Lyapunov exponent which may be 
understood by considering several probabilistic CA. We restrict our study to one dimensional, 
totalistic Boolean cellular automata with four, five and six inputs, since their number is 
reasonably manageable and their evolution can be efficiently implemented 

A Boolean CA F of range r is defined as a map on the set of configurations {x} with 
x = (xq, • • • , xn-i), Xi = 0, 1, and i = 0, . . . , N — 1 such that 



where x = x(t), x' = x(t + 1) and t = 0, 1 . . . . The map F is defined locally on every site i 



where {xi} r = (xj, . . . , Xj +r _i) is the neighborhood of range r of site i at time t, assuming 
periodic boundary conditions. For totalistic CA, the local function / is symmetric and 
depends only on s defined by 



x' = f(x), 

y' = (i-p)f(y)+pf(x), 



p c = 1 - exp(-A). 



(1) 



x' = F(x) 



by 



X'i = f({Xi}r) 



r-1 




j=0 
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That is x\ = f(s({xi} r )). It is useful to introduce the following operations between Boolean 
quantities: the sum modulo two (XOR), denoted by the symbol ©, and the AND operation, 
which is analogous to the usual multiplication and shares the same symbol. These operations 
can be performed between two configurations component by component. We introduce the 
difference, or damage, z = x © y, whose evolution is given by z' = F(x) © F(y) and we 
define the norm of z as |z| = (l/N) ^ Xi © yi. 

A function f(xi,...,Xj,..., Xi +r ) is sensitive to its j-th argument for a given neighbor- 
hood {{xi} r ) if the Boolean derivative 



df 



— f[ x ij ■ ■ ■ ) x ji ■ ■ ■ ) © f\ x i-> ■ ■ ■ j x j © 1) ■ ■ ■ ) 

{Xi}r 



dxj 

is 1. The Jacobian matrix J of F is an N x N matrix with components 

df 



J - (X) dx 3 



The matrix J is circular with zeroes everywhere except possibly on the main diagonal and 
the following r — 1 upper diagonals. 

It is possible to "Taylor expand" a Boolean function around a given point using Boolean 
derivatives [|IIJ]. To first order in |z| we have 

F(y) = F(x) © J(x) z (2) 

where denotes the Boolean multiplication of a matrix by a vector. Compared to algebraic 
multiplication of a matrix by a vector, the sum and multiplication of scalars are replaced 
by the XOR and the AND operations respectively. Using Eq. (0) we may approximate the 
evolution of the damage configuration z by 

z' = J(x) z. 

However, |z| grows at most linearly with t since a damage cannot spread to more than r 
neighbors in one time step: a fixed site i at time t + 1 can be damaged if at least one of 
its r neighbors at time t is damaged, but if more than one of the neighbors is damaged, the 
damage may cancel. Since 

i+r— 1 



%i — t I 3 Ji.n[X-)Zji 



z[ = 1 if Jjj (x)^ = 1 on an odd number of sites. In order to account for all possible damage 
spreading we choose to consider each damage independently. If, at time t, m damaged sites 
are present, we consider m replicas each one with a different damaged site. On each replica, 
the damage evolves for one time step, without interference effects and so on. 

This procedure is equivalent to choosing a vector £(0) = z(0) that evolves in time 
according to 

r = j(x)€ (3) 
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where now the matrix multiplication is algebraic. The components £j are positive integers 
that count the number of ways in which the initial damage can spread to site i at time t on 
the ensemble of replicas. 

We define the maximum Lyapunov exponent A of the cellular automaton F by 

1, 



A(x ) = lim lim — lo ft , 

where |£| may be taken as the Euclidean norm or as the sum of its components. The 
geometrical average /i of ones in the Jacobian matrix J is defined by 

T-l 



^(x ) = lim lim ( TT - V J, 




\ t=0 i,j 

We show in Fig. [I] the points (r/i, A) of all the totalistic CA with r = 4,5, and 6 that show 
complex non periodic space time patterns. 

The process defined by Eq. (|3|) may be viewed as a deterministic directed bond perco- 
lation problem where a site i at depth t is wet if > 0. The bonds exist where the 
components of J are 1. A first approximation is obtained by replacing J with a random 
matrix whose elements are zero except on the diagonal and the r — 1 following upper diago- 
nals, where they are equal to one with probability /i f] ; ll,12]. There is a critical value // c (r) 



below which the bond percolation process falls into the absorbing state so that the maximum 
eigenvalue of the product of random matrices is zero. We can further introduce a mean field 
approximation to the directed bond percolation process, which exhibits discrepancies only 
very near to \i c . In this case one can show that 

A = ln(r/x) (4) 

is an upper bound to the MLE of the product of random matrices. The behavior of A is 
plotted in Fig. [I], and we note that it is a good approximation to the generic behavior of CA. 
We also report the value r/j, c (r) for which the maximum eigenvalue of the product of random 
matrices is zero, corresponding to the percolation threshold for the directed bond percolation 
problem. We found that r/x c (r) ~ 1.3 regardless of r, a fact that can be understood from 
the following argument. The percolation cluster has on the average k bonds per site with 
k(t) = r/i c (r); for this percolation model, one connection to a wet neighbor at the previous 
time step is sufficient, and further connections do not alter wetting. 

We now discuss a synchronization mechanism for CA. Starting with two initial configu- 
rations chosen at random x(0) and y(0) we propose that 

x' = F(xV 



y' = S t (p)F(y)®5 t (p)F(x) l 

where S t (p) is a Boolean random diagonal matrix with elements sj(p) that, at each time 
step, take the value one with probability p and zero with probability 1 — p\ S(p) = I — S(p) 
and / is the identity matrix. On the average, y[ will be set to the value of x\ = /({xj}) on a 
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fraction p of sites. In this way we introduced a stochastic synchronization mechanism over 
a deterministic process. This stochatic mechanism can be considered as a "random field" 
approximation of an intermittent coupling generated by a deterministic chaotic process. The 
evolution equation for the difference z = x © y is 



i! = S(p) [F(x) © F(y)] . (5) 

The control and order parameters are p and h(p) = lim^oo limA^oo |z(£)| respectively. We 
say that x, the driver, and y, the driven system, synchronize when h(p) = 0. For p = 
both systems evolve independently, while for p — 1 they synchronize in just one step; we 
expect then to find a synchronization threshold p c . This behavior is shared by all the CA 
with complex non periodic space time patterns. All others synchronize for p ~ 0. This can 
be conversely expressed by saying that all CA that synchronize with a nontrivial p c exhibit 
complex non periodic space time patterns. 

For totalistic linear rules, whose evolution is given by 



3=0 



the synchronization equation fl5|) is equivalent to the dilution (with probability 1 — p) of 
the rule. For r = 2 the dilution problem is equivalent to the line z = of Ref. fI5| , whose 
transition belongs to the universality class of directed percolation (DP) fCIl . For a similar 
synchronization mechanism, it has been recently claimed that the elementary CA rule 18 does 
not Q and does || belong to the DP universality class. The presence of a single absorbing 
state and the absence of other conserved quantities (i.e. number of kinks) strongly suggests 
that the synchronization transition belongs to the DP universality class |3|||. The goal of 
this work was not that of computing critical exponents. However, since the critical point 
was located by means of the scaling law for the density of defects (in order to minimize finite 
size and time effects), we had the opportunity of computing the magnetic exponent /3 for all 
the CA studied. Due to the large number of CA examined, this computation was performed 
semi-automatically, and the precision of the resulting exponent is quite low, nontheless 
we have not found any example of non-DP behavior. We performed bulk simulations for 
N = 2000 and T = 4000 and checked that the results do not change for doubling of halving 
these figures. 

In Fig. |2] we plot the points (p c , A) for all totalistic CA with r = 4,5,6 and a non trivial 
value of p c . All the numerical experiments were performed using a parallel algorithm that 



takes care of all the values of p simultaneously |15 . 

Let us study the relation between p c and A by some random approximations. Near the 
synchronization threshold p c , we may expand y around x with the help of Eq. (0) so that 



z = S(p) J(x)z. (6) 

Eq. (0) may be written as 



z[ = Si(p) [Jj,i(x)^ © ■ ■ ■ © J M+r _i(x)2: i+r _i] . (7) 
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During the time evolution of a particular CA a fixed value of p is attained so that on average 
rp derivatives inside the parenthesis on the r.h.s. of Eq. (0) are different from zero. A first 
approximation, model A, is obtained by replacing the derivatives with random variables 
rrii(n) that are one with probability p and zero with probability (1 — p). That is 



A = Si(p) [mil^Zi © • • ■ © m i+r ^i(p)z i+r ^ 1 } . 

On every site and every time step, rp variables are chosen on average and if their sum is 
odd, then with probability p, z[ = 1. We then look for the synchronization threshold p c {p) 
and plot it as a function of A = ln(r/x), Eq. ([|). This is the curve labeled A in Fig. |2|. The 
predicted values of p c are generally higher than those found for CA for the same value of A 
possibly due to correlations among the derivatives. 

Since the typical "complex" CA pattern exhibit transient correlations ("triangles"), one 
can model them by choosing a fixed number k < r of derivatives equal to one. The simplest 
way is to take k = rp with rp an integer, model B, which . Then 



Z i = Si(p) [Zi © • • • © Zi+fc-l] • 

which is a dilution of the XOR with k inputs. This process is expected to belong to the 
same universality class of directed percolation |16]. The curve labeled B of Fig. |2] passes 



through all the calculated values for k = 2, ... 6. We note that this second model is a better 
fit of CA behavior. 

We can extend this last model allowing noninteger values of rp by 



Zi(t + 1) = s\{p)m\{p) © • • • © iZj+fc-i] 

which we call model In this way the average number of non zero derivatives is kp with 
< p < 1. Now A; is a free parameter, and this model can be useful to delimit the expected 
spread of (A,p c ) points. 

Since s\(p) and m\(p) are independent random variables we may write 



Zi(t + 1) = s\(q) © • • • © z i+ k-i] ■ 

where (1 — q) = (1 — p)p. In this guise, this is model B with k inputs. The synchronization 
threshold is given by p c (k, p) = 1 — (1 — q c (k))/p where q c {k) is the percolation threshold of 
the dilution of the XOR with k inputs. Using the approximation A = ln(kp) one has 

p c (M) = l-fc(l-g c (AO)exp(-A). (8) 

which bears some likeness to Eq. (|l|). The curves labeled Ci and C@ of Fig. |2] correspond 
to this last expression for k = 2 and k = 6 respectively. We note that the points (\,p c ) of 
almost all the CA considered fall between these two curves. 

For the totalistic one dimensional CA with r = 4, 5, 6 we can safely say that all CA 
with a positive p c exhibit complex non periodic space time patterns and conversely. These 
CA also have a positive MLE . We also showed that the synchronization of CA is a critical 
phenomenon similar to directed percolation. 

We proposed several approximations based on a combination of "linearization" of CA 
rules using Boolean Taylor expansions and stochasticity and showed the relation between the 
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synchronization threshold and the MLE. In particular, model C implies a relation similar 
to that found for continuous maps with the addition of a percolation threshold constant. 

An analogous mechanism can be applied to coupled map lattices; in this case p is the 
probability that a map Ui{t) takes the value Xi(t). One observes a synchronization transitions, 
but the critical value p c is not correlated to the usual MLE |HJ . 

We acknowledge partial financial support from project DGAPA-UNAM IN103595. 
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FIG. 1. The MLE of totalistic CA with 
r = 4, 5, 6 versus rfi. The continuous line repre- 
sents the mean field approximation A = ln(r/x). 
The dashed line marks the threshold rfi c . 
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FIG. 2. Relationship between p c and A for all 
CA with range r = 4, 5, 6 (markers) and complex 
space time patterns. The curves correspond to the 
various approximations, as specified in the legend. 
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